
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cumfreq, skew
from scipy.stats import norm

# 抽樣次數
N = 50000
# 樣本數
n_val = [10, 20, 30, 50, 100, 300, 500, 1000, 2000]
G1 = np.zeros((N, len(n_val)))

# 繪圖
plt.style.use('ggplot')
fig, axs = plt.subplots(9, 2, figsize = (12, 27))

for i, n in enumerate(n_val):
    # 生成 N rows n columns 服從標準常態分佈的樣本 for skewness
    x = norm.rvs(0, 1, (N, n))
    # 生成 G1 统计量
    G1 = np.sqrt(n / 6) * skew(x, axis = 1, bias = False) # compute skewness for each row

    # normal PDF
    x_pdf = np.linspace(min(G1), max(G1), 1000)
    y_pdf = norm.pdf(x_pdf, 0, 1)

    # 計算 G1 的 ECDF
    res = cumfreq(G1, numbins = N)
    ecdf_x = res.lowerlimit + np.linspace(0, res.binsize * res.cumcount.size, res.cumcount.size)
    cumcountprob = res.cumcount / N  # 累積次數除以總樣本數量 = 累積機率

    # normal CDF
    x_cdf = np.linspace(ecdf_x.min(), ecdf_x.max(), 1000)
    y_cdf = norm.cdf(x_cdf, 0, 1)

    # 設定子圖
    ax_pdf = axs[i, 0]
    ax_cdf = axs[i, 1]

    # G1 统计量的直方图與 normal PDF
    his = ax_pdf.hist(G1, bins = 100, density = True, color = '#99BBFF', edgecolor = '#FFFFFF'
                , label = 'G1 Statistic')
    ax_pdf.plot(x_pdf, y_pdf, color = '#FFA488', lw = 3, linestyle = '--', label = 'Z Distribution'
                , alpha = 0.9)
    ax_pdf.set_title('G1 Statistic: n = {}'.format(n), fontsize = 12)
    ax_pdf.set_ylabel('Density', fontsize = 12, color = 'black')
    legend0 = ax_pdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend0.get_frame().set_linestyle('-.')
    legend0.get_frame().set_alpha(0.4)
    ax_pdf.grid(True, linestyle = '-.')
    ax_pdf.tick_params(axis = 'both', labelsize = 11, colors = 'black')
    counts = his[0]
    ax_pdf.set_xlim([x_pdf.min() - 0.5, x_pdf.max() + 0.5])
    ax_pdf.set_ylim(-0.15, counts.max() + 0.05)

    # G1 统计量的盒鬚圖
    boxprops = dict(linestyle = '--', linewidth = 2, color = '#9999FF')
    whiskerprops = dict(color = '#9999FF', linewidth = 1.7, linestyle = '--')
    capprops = dict(color = '#9999FF', linewidth = 1.7)
    flierprops = dict(marker = 'o', markerfacecolor = '#FFFFBB', markersize = 4
                    , linestyle = 'none', markeredgecolor = 'none')   # define outliers
    medianprops = dict(color = '#FFFFBB', linewidth = 2, linestyle = '-')  # 设置中位数线的样式
    ax_box = ax_pdf.twinx()  # 使用 twinx 創建共享 x 軸的第二個 y 軸
    ax_box.boxplot(G1, vert = False, notch = True, boxprops = boxprops, flierprops = flierprops, whiskerprops = whiskerprops
                   , capprops = capprops, medianprops = medianprops, positions = [-0.12], widths = 0.03)
    ax_box.set_yticks([])  # 移除 y 軸刻度
    ax_box.set_ylim(-0.15, 0.15)  # 設置 y 軸範圍以顯示盒鬚圖

    # G1 统计量的 ECDF 與 normal CDF
    ax_cdf.plot(ecdf_x, cumcountprob, drawstyle = 'steps-pre', label = 'Empirical CDF'
                , color = '#AAAAAA', lw = 3)
    ax_cdf.plot(x_cdf, y_cdf, color = '#FFA488', label = 'Theoretical CDF', alpha = 0.9, lw = 3
                , linestyle = '--')
    ax_cdf.set_title('The Empirical CDF with {} steps'.format(n), fontsize = 12)
    ax_cdf.set_ylabel('Cumulative Probability', fontsize = 12, color = 'black')
    legend1 = ax_cdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend1.get_frame().set_linestyle('-.')
    legend1.get_frame().set_alpha(0.4)
    ax_cdf.grid(True, linestyle = '-.')
    ax_cdf.tick_params(axis = 'both', labelsize = 11, colors = 'black')

plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cumfreq, kurtosis
from scipy.stats import norm

# 抽樣次數
N = 50000
# 樣本數
n_val = [10, 20, 30, 50, 100, 300, 500, 1000, 3000, 5000]
G2 = np.zeros((N, len(n_val)))

# 繪圖
plt.style.use('ggplot')
fig, axs = plt.subplots(10, 2, figsize = (12, 30))

for i, n in enumerate(n_val):
    # 生成 N rows n columns 服從標準常態分佈的樣本 for kurtosis
    x = norm.rvs(0, 1, (N, n))
    # 生成 G2 统计量
    # G2 = np.sqrt(n / 24) * (kurtosis(x, axis = 1, bias = False, fisher = False) - 3)
    G2 = np.sqrt(n / 24) * (kurtosis(x, axis = 1, bias = False, fisher = True))

    # normal PDF
    x_pdf = np.linspace(min(G2), max(G2), 1000)
    y_pdf = norm.pdf(x_pdf, 0, 1)

    # 計算 G2 的 ECDF
    res = cumfreq(G2, numbins = N)
    ecdf_x = res.lowerlimit + np.linspace(0, res.binsize * res.cumcount.size, res.cumcount.size)
    cumcountprob = res.cumcount / N  # 累積次數除以總樣本數量 = 累積機率

    # normal CDF
    x_cdf = np.linspace(ecdf_x.min(), ecdf_x.max(), 1000)
    y_cdf = norm.cdf(x_cdf, 0, 1)

    # 設定子圖
    ax_pdf = axs[i, 0]
    ax_cdf = axs[i, 1]

    # G2 统计量的直方图與 normal PDF
    his = ax_pdf.hist(G2, bins = 100, density = True, color = '#F0ABF5', edgecolor = '#FFFFFF'
                , label = 'G2 Statistic')
    ax_pdf.plot(x_pdf, y_pdf, color = '#FFA488', lw = 3, linestyle = '--', label = 'Z Distribution'
                , alpha = 0.9)
    ax_pdf.set_title('G2 Statistic: n = {}'.format(n), fontsize = 12)
    ax_pdf.set_ylabel('Density', fontsize = 12, color = 'black')
    legend0 = ax_pdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend0.get_frame().set_linestyle('-.')
    legend0.get_frame().set_alpha(0.4)
    ax_pdf.grid(True, linestyle = '-.')
    counts = his[0]
    ax_pdf.set_xlim([x_pdf.min() - 0.5, x_pdf.max() + 0.5])
    ax_pdf.set_ylim(-0.15, counts.max() + 0.05)

    # G2 统计量的盒鬚圖
    boxprops = dict(linestyle = '--', linewidth = 2, color = '#9999FF')
    whiskerprops = dict(color = '#9999FF', linewidth = 1.7, linestyle = '--')
    capprops = dict(color = '#9999FF', linewidth = 1.7)
    flierprops = dict(marker = 'o', markerfacecolor = '#FFFFBB', markersize = 4
                    , linestyle = 'none', markeredgecolor = 'none')   # define outliers
    medianprops = dict(color = '#FFFFBB', linewidth = 2, linestyle = '-')  # 设置中位数线的样式
    ax_box = ax_pdf.twinx()  # 使用 twinx 創建共享 x 軸的第二個 y 軸
    ax_box.boxplot(G2, vert = False, notch = True, boxprops = boxprops, flierprops = flierprops, whiskerprops = whiskerprops
                   , capprops = capprops, medianprops = medianprops, positions = [-0.12], widths = 0.03)
    ax_box.set_yticks([])  # 移除 y 軸刻度
    ax_box.set_ylim(-0.15, 0.15)  # 設置 y 軸範圍以顯示盒鬚圖

    # G2 统计量的 ECDF 與 normal CDF
    ax_cdf.plot(ecdf_x, cumcountprob, drawstyle = 'steps-pre', label = 'Empirical CDF'
                , color = '#AAAAAA', lw = 3)
    ax_cdf.plot(x_cdf, y_cdf, color = '#FFA488', label = 'Theoretical CDF', alpha = 0.9, lw = 3
                , linestyle = '--')
    ax_cdf.set_title('The Empirical CDF with {} steps'.format(n), fontsize = 12)
    ax_cdf.set_ylabel('Cumulative Probability', fontsize = 12, color = 'black')
    legend1 = ax_cdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend1.get_frame().set_linestyle('-.')
    legend1.get_frame().set_alpha(0.4)
    ax_cdf.grid(True, linestyle = '-.')

plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cumfreq, skew, kurtosis
from scipy.stats import norm, chi2

# 抽樣次數
N = 50000
# 樣本數
n_val = [10, 20, 30, 50, 100, 300, 500, 1000, 2000]
G3 = np.zeros((N, len(n_val)))

# 繪圖
plt.style.use('ggplot')
fig, axs = plt.subplots(9, 2, figsize = (12, 27))

for i, n in enumerate(n_val):
    # 生成 N rows n columns 服從標準常態分佈的樣本 for kurtosis
    x = norm.rvs(0, 1, (N, n))
    # 生成 G3 统计量
    s = skew(x, axis = 1, bias = False)
    k = kurtosis(x, axis = 1, bias = False, fisher = True)
    G3 = (n / 6) * ((s ** 2) + ((k ** 2) / 4))

    # chi(2) PDF
    x_pdf = np.linspace(min(G3), max(G3), 1000)
    y_pdf = chi2.pdf(x_pdf, 2)

    # 計算 G3 的 ECDF
    res = cumfreq(G3, numbins = N)
    ecdf_x = res.lowerlimit + np.linspace(0, res.binsize * res.cumcount.size, res.cumcount.size)
    cumcountprob = res.cumcount / N  # 累積次數除以總樣本數量 = 累積機率

    # chi(2) CDF
    x_cdf = np.linspace(ecdf_x.min(), ecdf_x.max(), 1000)
    y_cdf = chi2.cdf(x_cdf, 2)

    # 設定子圖
    ax_pdf = axs[i, 0]
    ax_cdf = axs[i, 1]

    # G3 统计量的直方图與 chi(2) PDF
    his = ax_pdf.hist(G3, bins = 150, density = True, color = '#99d8c9', edgecolor = '#FFFFFF'
                , label = 'G3 Statistic')
    ax_pdf.plot(x_pdf, y_pdf, color = '#FFA488', lw = 3, linestyle = '--', label = '$\chi^2(2)$ Distribution'
                , alpha = 0.9)
    ax_pdf.set_title('G3 Statistic: n = {}'.format(n), fontsize = 12)
    ax_pdf.set_ylabel('Density', fontsize = 12, color = 'black')
    legend0 = ax_pdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend0.get_frame().set_linestyle('-.')
    legend0.get_frame().set_alpha(0.4)
    ax_pdf.grid(True, linestyle = '-.')
    counts = his[0]
    ax_pdf.set_xlim([-0.1, 11])
    ax_pdf.set_ylim(-0.15, counts.max() + 0.1)

    # G3 统计量的盒鬚圖
    boxprops = dict(linestyle = '--', linewidth = 2, color = '#9999FF')
    whiskerprops = dict(color = '#9999FF', linewidth = 1.7, linestyle = '--')
    capprops = dict(color = '#9999FF', linewidth = 1.7)
    flierprops = dict(marker = 'o', markerfacecolor = '#FFFFBB', markersize = 4
                    , linestyle = 'none', markeredgecolor = 'none')   # define outliers
    medianprops = dict(color = '#FFFFBB', linewidth = 2, linestyle = '-')  # 设置中位数线的样式
    ax_box = ax_pdf.twinx()  # 使用 twinx 創建共享 x 軸的第二個 y 軸
    ax_box.boxplot(G3, vert = False, notch = True, boxprops = boxprops, flierprops = flierprops, whiskerprops = whiskerprops
                   , capprops = capprops, medianprops = medianprops, positions = [-0.12], widths = 0.03)
    ax_box.set_yticks([])  # 移除 y 軸刻度
    ax_box.set_ylim(-0.15, 0.15)  # 設置 y 軸範圍以顯示盒鬚圖

    # G3 统计量的 ECDF 與 chi(2) CDF
    ax_cdf.plot(ecdf_x, cumcountprob, drawstyle = 'steps-pre', label = 'Empirical CDF'
                , color = '#AAAAAA', lw = 3)
    ax_cdf.plot(x_cdf, y_cdf, color = '#FFA488', label = 'Theoretical CDF', alpha = 0.9, lw = 3
                , linestyle = '--')
    ax_cdf.set_xlim([0, 12])
    ax_cdf.set_title('The Empirical CDF with {} steps'.format(n), fontsize = 12)
    ax_cdf.set_ylabel('Cumulative Probability', fontsize = 12, color = 'black')
    legend1 = ax_cdf.legend(edgecolor = '#666666', prop = {'size': 11})
    legend1.get_frame().set_linestyle('-.')
    legend1.get_frame().set_alpha(0.4)
    ax_cdf.grid(True, linestyle = '-.')

plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import skew, kurtosis, norm, chi2, uniform, t, f

# 1. 計算檢定統計量 G3 的函數 G3_stats(X)
def G3_stats(X):
    """
    Input:
        X: A numpy array of shape (N, n) where N is the number of samples and n is the size of each sample
    Output:
        G3: A numpy array of shape (N,) that contains the G3 statistic for each sample
        P_Val: A numpy array of shape (N,) that contains the p-value for each sample
    """
    n = X.shape[1]
    G1 = np.sqrt(n / 6) * skew(X, axis = 1, bias = False) # bias = False 表示 Unbiased
    G2 = np.sqrt(n / 24) * (kurtosis(X, axis = 1, bias = False, fisher = True))
    G3 = G1**2 + G2**2
    P_Val = 1 - chi2.cdf(G3, df = 2) # chi2.cdf 是指從 0 到 G3 的面積，1 - chi2.cdf 是指從 G3 到無限大的面積
    return G3, P_Val

distributions = {
    'N(0, 1)': lambda size: norm.rvs(loc = 0, scale = 1, size = size),
    'T(3)': lambda size: t.rvs(df = 3, size = size),
    'T(10)': lambda size: t.rvs(df = 10, size = size),
    'T(30)': lambda size: t.rvs(df = 30, size = size),
    'Unif(0, 1)': lambda size: uniform.rvs(loc = 0, scale = 1, size = size),
    '$\chi^2(4)$': lambda size: chi2.rvs(df = 4, size = size),
    '$\chi^2(8)$': lambda size: chi2.rvs(df = 8, size = size),
    '$\chi^2(20)$': lambda size: chi2.rvs(df = 20, size = size),
    'F(10, 10)': lambda size: f.rvs(dfn = 10, dfd = 10, size = size),
    'F(50, 50)': lambda size: f.rvs(dfn = 50, dfd = 50, size = size)
}


# 2. 快速選擇所需分配的函數 dist(dist_name, N, n)
def dist(dist_name, N, n):
    """
    Input:
        dist_name: A string that represents the distribution name
        N: The number of samples
        n: The size of each sample
    Output:
        distributions[dist_name]((N, n)): A numpy array of shape (N, n) that contains the random samples
    """
    return distributions[dist_name]((N, n))


# 設定參數
N = 50000 # number of samples
n = [10, 20, 30, 50, 100, 300, 500, 1000] # sample size
alpha = 0.05

# 繪圖設定
plt.style.use('ggplot')
fig, axs = plt.subplots(10, 2, figsize = (12, 30))
color = sns.color_palette('tab10', 10)
boxcolor = sns.color_palette('Set3', 8)
# color.pop(8)  # 移除第 9 个颜色，因为 F(50, 50) 的颜色太深
ylim_lower = [0.04, 0, 0, 0, -0.05, 0, 0, 0, 0, 0]
ylim_upper = [0.06, 1.15, 1.15, 0.35, 1.15, 1.15, 1.15, 1.15, 1.15, 1.15]

# 開始繪圖
for j, dist_name in enumerate(distributions.keys()):
    # 趨勢圖
    power = np.zeros(len(n))
    p_vals = [[] for _ in range(len(n))]  # 创建一个包含 len(n) 个空列表的列表
    for i in range(len(n)):
        x = dist(dist_name, N, n[i])
        G3, P_Val = G3_stats(x)
        power[i] = np.mean(P_Val < alpha)
        p_vals[i].extend(P_Val)

    axs[j, 0].plot(power, marker = 's', linestyle = '-', color = color[j], lw = 2.5, label = dist_name)
    axs[j, 0].axhline(y = alpha, color = '#FF7744', linestyle = '--', lw = 2, alpha = 0.7
                   , label = '$\\alpha = {}$'.format(alpha))
    axs[j, 0].set_xticks(np.arange(len(n)))
    axs[j, 0].set_xticklabels(n)
    axs[j, 0].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[j, 0].set_ylabel('Power', fontsize = 12, color = 'black')
    axs[j, 0].set_title('Population: ' + dist_name + ' at $\\alpha = {}$'.format(alpha), fontsize = 14, color = 'black')
    axs[j, 0].grid(True, linestyle = '-.')
    axs[j, 0].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    axs[j, 0].set_ylim([ylim_lower[j], ylim_upper[j]])
    legend = axs[j, 0].legend(prop = {'size': 10}, frameon = True, edgecolor = '#666666', loc = 'upper left')
    legend.get_frame().set_linestyle('-.')
    legend.get_frame().set_linewidth(1)
    legend.get_frame().set_alpha(0.8)

    # 将每个样本大小的 power_val 转换为 numpy 数组
    p_vals = [np.array(pv) for pv in p_vals]
    # 盒鬚圖
    for k in range(len(n)):
        boxprops = dict(linestyle = '-', linewidth = 1.7, color = 'k', facecolor = boxcolor[k])
        whiskerprops = dict(color = 'k', linewidth = 1.3, linestyle = '--')
        capprops = dict(color = 'k', linewidth = 1.3)
        flierprops = dict(marker = 'o', markerfacecolor = boxcolor[k], markersize = 6, linestyle = 'none'
                          , markeredgecolor = 'none')  # define outliers
        medianprops = dict(color = 'k', linewidth = 1.3, linestyle = '-')  # 设置中位数线的样式
        axs[j, 1].boxplot(p_vals[k], vert = True, notch = False, patch_artist = True, boxprops = boxprops
                          , flierprops = flierprops, whiskerprops = whiskerprops,capprops = capprops
                          , medianprops = medianprops, positions = [k], widths = 0.6)  # 使用单一位置值
    axs[j, 1].set_xticks(np.arange(len(n)))
    axs[j, 1].set_xticklabels(n)
    axs[j, 1].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[j, 1].set_ylabel('P Value', fontsize = 12, color = 'black')
    axs[j, 1].set_title('Boxplot of P-Values for ' + dist_name + ' at $\\alpha = {}$'.format(alpha), fontsize = 14, color = 'black')
    axs[j, 1].grid(True, linestyle = '-.')
    axs[j, 1].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    

plt.tight_layout()
plt.show()


# 設定參數
N = 50000 # number of samples
n = [10, 20, 30, 50, 100, 300, 500, 1000] # sample size
alpha = 0.01

# 繪圖設定
plt.style.use('ggplot')
fig, axs = plt.subplots(10, 2, figsize = (12, 30))
color = sns.color_palette('tab10', 10)
boxcolor = sns.color_palette('Set3', 8)
ylim_lower = [0.04, 0, 0, 0, -0.05, 0, 0, 0, 0, 0]
ylim_upper = [0.06, 1.15, 1.15, 0.35, 1.15, 1.15, 1.15, 1.15, 1.15, 1.15]

# 開始繪圖
for j, dist_name in enumerate(distributions.keys()):
    # 趨勢圖
    power = np.zeros(len(n))
    p_vals = [[] for _ in range(len(n))]  # 创建一个包含 len(n) 个空列表的列表
    for i in range(len(n)):
        x = dist(dist_name, N, n[i])
        G3, P_Val = G3_stats(x)
        power[i] = np.mean(P_Val < alpha)
        p_vals[i].extend(P_Val)

    axs[j, 0].plot(power, marker = 's', linestyle = '-', color = color[j], lw = 2.5, label = dist_name)
    axs[j, 0].axhline(y = alpha, color = '#FF7744', linestyle = '--', lw = 2, alpha = 0.7
                   , label = '$\\alpha = {}$'.format(alpha))
    axs[j, 0].set_xticks(np.arange(len(n)))
    axs[j, 0].set_xticklabels(n)
    axs[j, 0].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[j, 0].set_ylabel('Power', fontsize = 12, color = 'black')
    axs[j, 0].set_title('Population: ' + dist_name + ' at $\\alpha = {}$'.format(alpha), fontsize = 14, color = 'black')
    axs[j, 0].grid(True, linestyle = '-.')
    axs[j, 0].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    axs[j, 0].set_ylim([ylim_lower[j], ylim_upper[j]])
    legend = axs[j, 0].legend(prop = {'size': 10}, frameon = True, edgecolor = '#666666', loc = 'upper left')
    legend.get_frame().set_linestyle('-.')
    legend.get_frame().set_linewidth(1)
    legend.get_frame().set_alpha(0.8)

    # 将每个样本大小的 power_val 转换为 numpy 数组
    p_vals = [np.array(pv) for pv in p_vals]
    # 盒鬚圖
    for k in range(len(n)):
        boxprops = dict(linestyle = '-', linewidth = 1.7, color = 'k', facecolor = boxcolor[k])
        whiskerprops = dict(color = 'k', linewidth = 1.3, linestyle = '--')
        capprops = dict(color = 'k', linewidth = 1.3)
        flierprops = dict(marker = 'o', markerfacecolor = boxcolor[k], markersize = 6, linestyle = 'none'
                          , markeredgecolor = 'none')  # define outliers
        medianprops = dict(color = 'k', linewidth = 1.3, linestyle = '-')  # 设置中位数线的样式
        axs[j, 1].boxplot(p_vals[k], vert = True, notch = False, patch_artist = True, boxprops = boxprops
                          , flierprops = flierprops, whiskerprops = whiskerprops,capprops = capprops
                          , medianprops = medianprops, positions = [k], widths = 0.6)  # 使用单一位置值
    axs[j, 1].set_xticks(np.arange(len(n)))
    axs[j, 1].set_xticklabels(n)
    axs[j, 1].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[j, 1].set_ylabel('P Value', fontsize = 12, color = 'black')
    axs[j, 1].set_title('Boxplot of P-Values for ' + dist_name + ' at $\\alpha = {}$'.format(alpha), fontsize = 14, color = 'black')
    axs[j, 1].grid(True, linestyle = '-.')
    axs[j, 1].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    

plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import skew, kurtosis, norm, chi2, uniform, t, f

# 1. 計算檢定統計量 G3 的函數 G3_stats(X)
def G3_stats(X):
    """
    Input:
        X: A numpy array of shape (N, n) where N is the number of samples and n is the size of each sample
    Output:
        G3: A numpy array of shape (N,) that contains the G3 statistic for each sample
        P_Val: A numpy array of shape (N,) that contains the p-value for each sample
    """
    n = X.shape[1]
    G1 = np.sqrt(n / 6) * skew(X, axis = 1, bias = False) # bias = False 表示 Unbiased
    G2 = np.sqrt(n / 24) * (kurtosis(X, axis = 1, bias = False, fisher = True))
    G3 = G1**2 + G2**2
    P_Val = 1 - chi2.cdf(G3, df = 2) # chi2.cdf 是指從 0 到 G3 的面積，1 - chi2.cdf 是指從 G3 到無限大的面積
    return G3, P_Val

distributions = {
    'N(0, 1)': lambda size: norm.rvs(loc = 0, scale = 1, size = size),
    'T(3)': lambda size: t.rvs(df = 3, size = size),
    'T(10)': lambda size: t.rvs(df = 10, size = size),
    'T(30)': lambda size: t.rvs(df = 30, size = size),
    'Unif(0, 1)': lambda size: uniform.rvs(loc = 0, scale = 1, size = size),
    '$\chi^2(4)$': lambda size: chi2.rvs(df = 4, size = size),
    '$\chi^2(8)$': lambda size: chi2.rvs(df = 8, size = size),
    '$\chi^2(20)$': lambda size: chi2.rvs(df = 20, size = size),
    'F(10, 10)': lambda size: f.rvs(dfn = 10, dfd = 10, size = size),
    'F(50, 50)': lambda size: f.rvs(dfn = 50, dfd = 50, size = size)
}


# 2. 快速選擇所需分配的函數 dist(dist_name, N, n)
def dist(dist_name, N, n):
    """
    Input:
        dist_name: A string that represents the distribution name
        N: The number of samples
        n: The size of each sample
    Output:
        distributions[dist_name]((N, n)): A numpy array of shape (N, n) that contains the random samples
    """
    return distributions[dist_name]((N, n))


import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import jarque_bera, normaltest, kstest, shapiro, anderson, cramervonmises
from statsmodels.stats.diagnostic import lilliefors
import warnings

# 忽略特定警告
warnings.filterwarnings("ignore")

# 設定參數
N = 50000  # number of samples
n = [10, 20, 30, 50, 100, 300, 500, 1000]  # sample size
alpha = 0.05

# 繪圖設定
plt.style.use('ggplot')
fig, axs = plt.subplots(10, 1, figsize = (12, 30))
def distplot(power, color, label):
    axs[k].plot(power, marker = 's', linestyle = '-', color = color, lw = 2.5, label = label)

# 開始繪圖
for k, dist_name in enumerate(distributions.keys()):

    # 初始化 power
    power_jb = np.zeros(len(n))
    power_G3 = np.zeros(len(n))
    power_nt = np.zeros(len(n))
    power_ks = np.zeros(len(n))
    power_sw = np.zeros(len(n))
    power_ad = np.zeros(len(n))
    power_lf = np.zeros(len(n))
    power_cv = np.zeros(len(n))

    # 計算 power for each sample size
    for i in range(len(n)):
        # 選擇分佈
        x = dist(dist_name, N, n[i])

        # 自己的 j-b test
        G3, P_Val = G3_stats(x)
        # scipy j-b test
        p_val_jb = jarque_bera(x, axis = 1)[1]
        # scipy normaltest
        p_val_nt = normaltest(x, axis = 1)[1]
        # scipy kstest
        p_val_ks = np.zeros(N)
        for j in range(N):
            p_val_ks[j] = kstest(x[j], 'norm')[1]
        # scipy shapiro test
        p_val_sw = np.zeros(N)
        for j in range(N):
            p_val_sw[j] = shapiro(x[j])[1]
        # anderson-darling test
        stat_ad = np.zeros(N)
        critical_val_ad = np.zeros(N)
        for j in range(N):
            result = anderson(x[j], dist = 'norm')
            stat_ad[j] = result.statistic
            critical_val_ad[j] = result.critical_values[2]  # 第三個臨界值，為 5% 顯著性水平
        # lilliefors test
        p_val_lf = np.zeros(N)
        for j in range(N):
            p_val_lf[j] = lilliefors(x[j])[1]
        # cramér-von mises test
        p_val_cv = np.zeros(N)
        for j in range(N):
            result_cv = cramervonmises(x[j], 'norm') # 設定為 'norm' 表示比較常態分佈
            p_val_cv[j] = result_cv.pvalue

        # 計算 power
        power_G3[i] = np.mean(P_Val < alpha)
        power_jb[i] = np.mean(p_val_jb < alpha)
        power_nt[i] = np.mean(p_val_nt < alpha)
        power_ks[i] = np.mean(p_val_ks < alpha)
        power_sw[i] = np.mean(p_val_sw < alpha)
        power_ad[i] = np.mean(stat_ad > critical_val_ad)  # 檢驗統計量 > 5% 的臨界值就 reject H0
        power_lf[i] = np.mean(p_val_lf < alpha)
        power_cv[i] = np.mean(p_val_cv < alpha)

    # 繪圖
    axs = axs.ravel()
    power = [power_G3, power_jb, power_nt, power_ks, power_sw, power_ad, power_lf, power_cv]
    color = sns.color_palette('Set2', 8)
    label = ['G3 Statistic', 'Jarque-Bera Test', 'Pearson Normal Test', 'Kolmogorov-Smirnov Test'
                , 'Shapiro-Wilk Test', 'Anderson-Darling Test', 'Lilliefors Test'
                , 'Cramér-von Mises Test']
    # 依序繪製各種檢定的 power
    for c in range(len(power)):
        distplot(power[c], color[c], label[c])
    axs[k].axhline(y = alpha, color = '#FF7744', linestyle = '--', lw = 2, alpha = 0.8)
    axs[k].set_xticks(np.arange(len(n)))
    axs[k].set_xticklabels(n)
    axs[k].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[k].set_ylabel('Power', fontsize = 12, color = 'black')
    axs[k].set_title('Population: ' + dist_name + ' at $\\alpha = {}$'.format(alpha)
                     , fontsize = 15, color = 'black')
    axs[k].grid(True, linestyle = '-.')
    axs[k].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    legend = axs[k].legend(prop = {'size': 10}, frameon = True, edgecolor = '#666666'
                          , loc = 'upper right', bbox_to_anchor = (1.29, 0.95))
    legend.get_frame().set_linestyle('-.')
    legend.get_frame().set_linewidth(1)
    legend.get_frame().set_alpha(0.8)
    # axs[k].set_ylim([0, 0.1])

plt.tight_layout()
plt.show()
# 跑了 194 m 45.7 s


import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import jarque_bera, normaltest, kstest, shapiro, anderson, cramervonmises
from statsmodels.stats.diagnostic import lilliefors
import warnings

# 忽略特定警告
warnings.filterwarnings("ignore")

# 設定參數
N = 50000  # number of samples
n = [10, 20, 30, 50, 100, 300, 500, 1000]  # sample size
alpha = 0.01

# 繪圖設定
plt.style.use('ggplot')
fig, axs = plt.subplots(10, 1, figsize = (12, 30))
def distplot(power, color, label):
    axs[k].plot(power, marker = 's', linestyle = '-', color = color, lw = 2.5, label = label)

# 開始繪圖
for k, dist_name in enumerate(distributions.keys()):

    # 初始化 power
    power_jb = np.zeros(len(n))
    power_G3 = np.zeros(len(n))
    power_nt = np.zeros(len(n))
    power_ks = np.zeros(len(n))
    power_sw = np.zeros(len(n))
    power_ad = np.zeros(len(n))
    power_lf = np.zeros(len(n))
    power_cv = np.zeros(len(n))

    # 計算 power for each sample size
    for i in range(len(n)):
        # 選擇分佈
        x = dist(dist_name, N, n[i])

        # 自己的 j-b test
        G3, P_Val = G3_stats(x)
        # scipy j-b test
        p_val_jb = jarque_bera(x, axis = 1)[1]
        # scipy normaltest
        p_val_nt = normaltest(x, axis = 1)[1]
        # scipy kstest
        p_val_ks = np.zeros(N)
        for j in range(N):
            p_val_ks[j] = kstest(x[j], 'norm')[1]
        # scipy shapiro test
        p_val_sw = np.zeros(N)
        for j in range(N):
            p_val_sw[j] = shapiro(x[j])[1]
        # anderson-darling test
        stat_ad = np.zeros(N)
        critical_val_ad = np.zeros(N)
        for j in range(N):
            result = anderson(x[j], dist = 'norm')
            stat_ad[j] = result.statistic
            critical_val_ad[j] = result.critical_values[4]  # 第 5 個臨界值，為 1% 顯著性水平
        # lilliefors test
        p_val_lf = np.zeros(N)
        for j in range(N):
            p_val_lf[j] = lilliefors(x[j])[1]
        # cramér-von mises test
        p_val_cv = np.zeros(N)
        for j in range(N):
            result_cv = cramervonmises(x[j], 'norm') # 設定為 'norm' 表示比較常態分佈
            p_val_cv[j] = result_cv.pvalue

        # 計算 power
        power_G3[i] = np.mean(P_Val < alpha)
        power_jb[i] = np.mean(p_val_jb < alpha)
        power_nt[i] = np.mean(p_val_nt < alpha)
        power_ks[i] = np.mean(p_val_ks < alpha)
        power_sw[i] = np.mean(p_val_sw < alpha)
        power_ad[i] = np.mean(stat_ad > critical_val_ad)  # 檢驗統計量 > 1% 的臨界值就 reject H0
        power_lf[i] = np.mean(p_val_lf < alpha)
        power_cv[i] = np.mean(p_val_cv < alpha)

    # 繪圖
    axs = axs.ravel()
    power = [power_G3, power_jb, power_nt, power_ks, power_sw, power_ad, power_lf, power_cv]
    color = sns.color_palette('Set2', 8)
    label = ['G3 Statistic', 'Jarque-Bera Test', 'Pearson Normal Test', 'Kolmogorov-Smirnov Test'
                , 'Shapiro-Wilk Test', 'Anderson-Darling Test', 'Lilliefors Test'
                , 'Cramér-von Mises Test']
    # 依序繪製各種檢定的 power
    for c in range(len(power)):
        distplot(power[c], color[c], label[c])
    axs[k].axhline(y = alpha, color = '#FF7744', linestyle = '--', lw = 2, alpha = 0.8)
    axs[k].set_xticks(np.arange(len(n)))
    axs[k].set_xticklabels(n)
    axs[k].set_xlabel('Sample size', fontsize = 12, color = 'black')
    axs[k].set_ylabel('Power', fontsize = 12, color = 'black')
    axs[k].set_title('Population: ' + dist_name + ' at $\\alpha = {}$'.format(alpha)
                     , fontsize = 15, color = 'black')
    axs[k].grid(True, linestyle = '-.')
    axs[k].tick_params(axis = 'both', colors = 'black', labelsize = 11)
    legend = axs[k].legend(prop = {'size': 10}, frameon = True, edgecolor = '#666666'
                          , loc = 'upper right', bbox_to_anchor = (1.29, 0.95))
    legend.get_frame().set_linestyle('-.')
    legend.get_frame().set_linewidth(1)
    legend.get_frame().set_alpha(0.8)
    # axs[k].set_ylim([0, 0.1])

plt.tight_layout()
plt.show()

# 跑了 268 m 9.5 s

第四次作品：以蒙地卡羅實驗驗證 J-B 檢定統計量的卡方分配與檢定力的表現

作品目標：

目標一：通過蒙地卡羅實驗來驗證以 Skewness 為基礎的 $G_1$ 統計量，以 Kurtosis 為基礎的 $G_2$ 統計量是否服從標準常態 N(0, 1)，以及由兩者結合形成的 $G_3$ 統計量是否服從卡方 $\chi^2(2)$。

注意事項與討論：

注意事項與討論：

注意事項與討論：

目標二：從不同的母體分配中進行抽樣，並用 $G_3$ 檢定統計量以蒙地卡羅模擬方式來進行檢定，驗證其在來自包含常態分佈下的各種分佈，以及在不同樣本數下檢定力（Power）的表現。

設定： (1) 計算檢定統計量 $G_3$ 的函數 G3_stats(X)，以及 (2) 快速選擇所需分配的函數 dist(dist_name, N, n)

注意事項與討論：

說明：與前一實驗設定相同，但將 Type I error 設定改為 0.01，將小於 0.01 的 P 值取平均，估計出檢定力 Power，觀察在不同分佈下樣本數對於 Power 的影響。

注意事項與討論：

目標三：本目標延續目標二的實驗設定，並加上其他各種著名的常態檢定方法相互進行比較，觀察它們之間在不同樣本數 n 以及不同分佈下檢定力 Power 的表現及差異。

設定： (1) 計算檢定統計量 $G_3$ 的函數 `G3_stats(X)`，以及 (2) 快速選擇所需分配的函數 `dist(dist_name, N, n)`

注意事項與討論：

說明：與前一實驗設定相同，但將 Type I error 設定改為 0.01，觀察這些常態檢定方法之間在不同樣本數 n 以及不同分佈下檢定力 Power 的表現。

注意事項與討論：

附錄：